package scu.maqiang.piezoelectricity;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3T4PE;
import scu.maqiang.mesh.Mesh3T4;
import scu.maqiang.numeric.DirectSSolver;
import scu.maqiang.numeric.SRMatrix;

public class Piezoelectricity3DT4 {

	public static void main(String[] args) {
		long begin = System.nanoTime();
		Mesh3T4 mesh = new Mesh3T4().cube(10, 10, 10);
		FES3T4PE fs = new FES3T4PE(mesh);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		double[] RHS = new double[fs.GetNdof()];
		double[] x = RHS.clone();
		double C11 = 14.9e10;
		double C33 = 14.3e10;
		double C12 = 9.8e10;
		double C13 = 9.8e10;
		double G12 = 2.2e10;
		double G13 = 2.4e10;
		double e15 = 11.7;
		double e31 = -5.4;
		double e33 = 13.5;
		double ep12 = 8.0e-9;
		double ep33 = 7.2e-9;
		double[][] matD = {{C11, C12, C13,   0,   0,   0,    0,    0, -e31},
		                   {C12, C11, C13,   0,   0,   0,    0,    0, -e31}, 
		                   {C13, C13, C33,   0,   0,   0,    0,    0, -e33},
		                   {  0,   0,   0, G12,   0,   0,    0, -e15,    0},
		                   {  0,   0,   0,   0, G12,   0, -e15,    0,    0},
		                   {  0,   0,   0,   0,   0, G13,    0,    0,    0},
		                   {  0,   0,   0,   0, e15,   0, ep12,    0,    0}, 
		                   {  0,   0,   0, e15,   0,   0,    0, ep12,    0},
		                   {e31, e31, e33,   0,   0,   0,    0,    0, ep33}};
		fs.assembleStiff(new double[][][] {matD}, BVPType.CONSTITUTE_MATRIX_COMMON, A);
	    fs.assembleSource(new double[] {1000.0, 0.0, 0.0, 0.0}, BVPType.COMMON, RHS);
	    fs.applyBC_MBN(A, RHS, new int[] {0, 1, 2, 3}, new double[] {0.0, 0.0, 0.0, 0.0}, 4);
	    //A.save("D:\\EclipseWorkspace2\\PyFem\\Numeric\\PEmatFile.txt");
	    //MVO.saveArray("D:\\EclipseWorkspace2\\PyFem\\Numeric\\PERHS.txt", RHS);

		DirectSSolver dsolver = new DirectSSolver(A);

//	    //DirectSSolver dsolver = new DirectSSolver(A, fs.getnBand());
//	    //dsolver.PLUSolve(RHS, x, 1);
//	    NewIterSSolver solver = new NewIterSSolver(A);
//	    solver.GMRES_SSOR(RHS, x, 1.75, 400, 2);
//	    //solver.PCGSSSOR(RHS, x, 1.75, 2);
//	    
//	    double[][] uvwphi = new double[4][mesh.getNv()];
//	    fs.extractComponent(x, uvwphi);
//	    mesh.toTecplot("PESolution3D.dat", uvwphi);
//	    System.out.println("SOlution: ");
//	    System.out.println("u min:" + MVO.min(uvwphi[0]) + "\tmax: " + MVO.max(uvwphi[0]));
//	    System.out.println("v min:" + MVO.min(uvwphi[1]) + "\tmax: " + MVO.max(uvwphi[1]));
//	    System.out.println("v min:" + MVO.min(uvwphi[2]) + "\tmax: " + MVO.max(uvwphi[2]));
//	    System.out.println("phi min:" + MVO.min(uvwphi[3]) + "\tmax: " + MVO.max(uvwphi[3]));
		long end = System.nanoTime();
		System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
	}

}
